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The ‘standard’ model of cosmology is founded on the basis that the expansion rate of the universe is 
accelerating at present — as was inferred originally from the Hubble diagram of Type la snpernovae. 
There exists now a much bigger database of supernovae so we can perform rigorous statistical tests 
to check whether these ‘standardisable candles’ indeed indicate cosmic acceleration. Taking account 
of the empirical procedure by which corrections are made to their absolute magnitudes to allow for 
the varying shape of the light curve and extinction by dust, we find, rather surprisingly, that the 
data are still quite consistent with a constant rate of expansion. 


I. INTRODUCTION 

In the late 1990’s, studies of Type la supernovae (SN 
la) showed that the expansion rate of the universe ap¬ 
pears to be accelerating as if dominated by a cosmolog¬ 
ical constantPl^. Since then supernova cosmology has 
developed rapidly as an important probe of ‘dark en¬ 
ergy’. Empirical corrections are made to reduce the scat¬ 
ter in the observed magnitudes by exploiting the observed 
(anti) correlation bet ween the peak luminosity and the 
light curve widtfP^^. Other such correlations have since 
been found e.g. with the host galaxy mas^and metallic- 
itjPl. Cosmological parameters are then fitted, along with 
the parameters det ermining the light curves, by simple 
minimisatiorpJSHil] This metho d has a number of pit- 
falls as has been emphasised earlieifi^^i^. 

With ever increasing precision and size of SN la 
datasets, it is important to also improve the statistical 
analysis of the data. To accomodate model comparison, 
previous worl^l^ has introduced likelihood maximisa¬ 
tion. In this work we present an improved maximum 
likelihood analysis, finding rather different results. 

II. SUPERNOVA COSMOLOGY 

There are several approaches to making SN la ‘stan- 
dardiseable candles’. The different philosophies lead to 
mildly different results but the overall picture seems 
consistent^. In this paper we adopt the transparent 
approac h of ‘ Spectral Adaptive Lightcurve Template 2’ 
(SALT2)) ^^b9 | .^vtierein the SN la are standardised by fit¬ 
ting their light curve to an empirical template, and the 
parameters of this fit are used in the cosmological analy¬ 
sis. Every SN la is assigned three parameters, one being 
m*g, the apparent magnitude at maximum (in the rest 
frame ‘S-band’), while the other two describe the light 
curve shape and colour corrections: xi and c. The dis¬ 
tance modulus is then taken to be: 

= m*B - M + axi - j3c, (I) 

where M is the absolute magnitude, and a and /3 are 
assumed to be constants for all SN la. These global 


constants are fitted along with the cosmological param¬ 
eters. The physical mechanism(s) which give rise to the 
cor relatio ns that underlie these corrections remain uncer- 
tairPSHU. The SN la distance modulus is then compared 
to the expectation in the standard ACDM cosmological 
model: 

/r = 25 -I- 5 log 2 o(dL/Mpc), where: 

dii = c/Ho, i7o = lOO/i kms“^Mpc“^, 

H = i?0'\/^m(l + z)^ + -b (2) 

where di^,dB_,H are the luminosity distance, Hub¬ 
ble distance and Hubble parameter respectively, and 
are the matter, cosmological constant and 
curvature density in units of the critical densitj^^. There 
is a degeneracy between Hq and Mq so we fix the value 
of the Hubble parameter today to h = 0.7 which is con¬ 
sistent with independent measurements. 

III. MAXIMUM LIKELIHOOD ESTIMATORS 

To find the maximum likelihood estimator (MLE) from 
the data, we must define the appropriate likelihood: 

C = probability density(data|model), 

i.e. we have to first specify our model of the data. For 
a given SN la, the true data {m*g,Xi,c) are drawn from 
some global distribution. These values are contaminated 
by various sources of noise, yielding the observed values 
{m*g,xi,c). Assuming the SALT2 model is correct, only 
the true values obey equation Q. However when the 
experimental uncertainty is of the same order as the in¬ 
trinsic variance as in the present case, the observed value 
is not a good estimate of the true value. Parameterising 
the cosmological model by 6, the likelihood function can 
be written as: 

C= p[{fh*g,xi,c)\0] (3) 

p[{m*g, xi, c)\{M, xi, c), 6] p[{M, cci, c)|6l]dMda;idc, 
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FIG. 1. Distribution of the stretch and colour correction pa¬ 
rameters in the JLA sampl^^, with gaussians superimposed. 


which shows explicitly where the experimental uncertain¬ 
ties enter (first factor) and where the variances of the 
intrinsic distributions enter (second factor). 

Having a theoretically well-motivated distribution for 
the light curve parameters would be helpful, however this 
is not available. For simplicity we adopt global, indepen¬ 
dent gaussian distributions for all parameters, M, xi and 
c (see Fig. [^, i.e. model their probability density as: 

p[{M,xi,c)\9] = p{M\9)p{xi\9)p{c\9), where: 
p{M\9) = (27r(T^J"^/^ exp [(M - Mq) /(Tmo? /sj , 

p{xi\9) = (27rcr2^_J-^/2gxp|_ [(xi -a;i,o)/cr^i,o]V2} , 
p{c\9) = exp [(c - Cq) /(Jcaf /2} . (4) 

All 6 free parameters {Ma, ctmoi 2 : 1 , 0 , cq, CTco} are 
fitted along with the cosmological parameters and we 
include them in 9. Introducing the vectors Y = 
{Ml, xii. Cl,... Mjv, xiAT, Cat}, the zero-points Yb, and 
the matrix S; = diag^a]^^, g, a'^^, ■ ■ ■), the probability 
density of the true parameters writes: 

p{Y\9) = \2nYi\-^/\xp[-{Y-Yo)E-\Y-Yo f/2] , 

( 5 ) 


where | ... | denotes the determinant of a matrix. What 
remains is to specify the model of uncertainties on 
the data. Introducing another set of vectors X = 
{mg^,xii,ci,...}, the observed X, and the estimated 
experimental covariance matrix So (including both sta¬ 
tistical and systematic errors), the probability density of 
the data given some set of true parameters is: 


p{X\X,9) 


|27rEd| exp 


-iX-X)Y-\X-Xf/2 


To combine the exponentials we introduce the vector Z = 
{m|ji — ^ 1 , ill. Cl,... } and the block diagonal matrix 


A = 



V 


0 0 
1 0 
0 1 

0 


\ 

0 


( 7 ) 


With these, we have X — X = (ZA ^ — Y)A and so 
p{X\X,9) = p{Z\Y,9). The likelihood is then 


C= jp{Z\Y,9) p{Y\9)dY (8) 

= j dY 

X exp (-(y - Yo)YY\Y - Yof/2) 

X exp (-(F - ZA-^)AY-^A^{Y - ZA-^f , 


which can be integrated analytically to obtain: 


£= \2TT{Yd + A^YiA)\ 


-1/2 


( 9 ) 


X exp 


-{Z - YoA)iYd + A^YiA)-^{Z - Y^AY/2 


This is the likelihood (equation ([^) for the simple model 
of equation Q, and the quantity which we maximise 
in order to derive confidence limits. The 10 parameters 
we fit are (Hm, Ha, a, xpo, q,/3, cq, CTcq, Mq, ctmoI- We 
stress that it is necessary to consider all of these together 
and Hm and Ha have no special status in this regard. The 
advantage of our method is that we get a goodness-of-fit 
statistic in the likelihood which can be used to compare 
models or judge whether a particular model is a good 
fit. Note that the model is not just the cosmology, but 
includes modelling the distributions of xi and c. 

With this MLE, we can construct a confidence region 
in the 10-dimensional parameter space by defining its 
boundary as one of constant L. So long as we do not 
cross a boundary in parameter space, this volume will 
asymptotically have the coverage probability 


n-21og£/£„ 


Pcov — 


f^2{x-,v)dx, 


( 10 ) 


where /^2 [x] v) is the pdf of a chi-squared random vari¬ 
able with v degrees of freedom, and £max is the max¬ 
imum likelihood. With 10 parameters in the present 
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model, the values Pcov — {0.68 (“lcr”),0.95 (“2 (t”)} give 
—21og£/£max — {11.54,18.61} respectively. 

To eliminate the so-called ‘nuisance parameters’, we 
set similar bounds on the profile likelihood. Writing the 
interesting parameters as 9 and nuisance parameters as 
(j), the profile likelihood is defined as 

£p(0) = max£(0, 0). (11) 

<t> 


We substitute C by £p in equation (101 in order to 
construct confidence regions in this lower dimensional 
space; v is now the dimension of the remaining param¬ 
eter space. Looking at the flm — SIa plane, we have for 
Pcov ^ {0.68 (“lcr”),0.95 (“2 ct”),0.997 (“3ct”)}, the val¬ 
ues —21og£p/£max — {2.30,6.18,11.8} respectively. 



A. Comparison to other methods 

It is illuminating to relate our work to previously used 
methods in SN la analyses. One methocP^l maximises a 
likelihood, which is written in the case of uncorrelated 
magnitudes as 

C oc ]J(27rCT2^t)-i/2 exp > (12) 

so it integrates over ^sn to unity and can be used for 
model comparison. From Equation ([^ we see that this 
corresponds to assuming flat distributions for Xi and c. 
However the actual distributions of xi and c are close to 
gaussian, as seen in Fig.[^ Moreover although this likeli¬ 
hood apparently integrates to unity, it accounts for only 
the rrfg data. Integration over the Xi,c data demands 
compact support for the flat distributions so the normal¬ 
isation of the likelihood becomes arbitrary, making model 

comparison tricky. _ 

More commonly usedl^ is the ‘constrained 

= + o-hit), (13) 

but this cannot be used to compare models, since it 
is tuned to be 1 per degree of freedom for the ACDM 
model by adjusting an arbitrary er ror g jnt added to each 
data point. This has been criticisecP^E^, nevertheless the 
method continues to be widely used and the results pre¬ 
sented without emphasising that it is intended only for 
parameter estimation for the assumed (ACDM) model, 
rather than determining if this is indeed the best model. 

IV. ANALYSIS OF JLA CATALOGUE 


FIG. 2. Contour plot of the profile likelihood in the Dm — Da 
plane. We show 1, 2 and Su contonrs, regarding all other 
parameters as nnisance parameters, as red dashed lines, while 
the blue lines are 1 and 2a contours from the 10-dimensional 
parameter space projected on to this plane. 


c are well modelled as gaussians. Maximisation of the 
likelihood under specific constraints is summarised in 
Table |l] and the profile likelihood contours in the Dm —Da 
plane are shown in Fig. In Fig. we compare the 
measured distance modulus, /isN = — Mq -I- axi — (3c 

with its expected value in two cosmological models: 
‘ACDM’ is the best fit accelerating universe while ‘Milne’ 
is an universe expanding with constant velocity. The 
error bars are the square root of the diagonal elements 
of Ei -h AT-iEdA-i so include both experimental 
uncertainties and intrinsic dispersion. We show also the 
residuals with respect to the Milne model (which has 
been raised to take into account the change in Mq). 

To assess how well our model describes the data, we 
present in Fig.j^the ‘pull’ distribution. These are defined 
as the normalised, decorrelated residuals of the data, 

pulls = (Z - Yo^)C^”\ (14) 

where U is the upper triangular Cholesky factor of the 
covariance matrix Ej -|- A'^EjA. Performing a K-S test, 
comparing the pull distribution to a unit variance gaus¬ 
sian gives a p-value of 0.1389. 

To check the validity of our method and approxima¬ 
tions, we do a Monte Carlo simulation of experimental 
outcomes from a model with parameters matching our 
best fit (see Table Figure [5] shows the distribution of 
—21og[£true/'Cmax]j which is just as is expected. 


We focus on the Joint Lightcurve Analysis (JLA) 
catalogu^A^. (All data used are available on http: // 
supernovae.in2p3.fr/sdss_snls_jla/ReadMe.html 
— we use the covmat_v6.) As shown already in Fig. 
the distributions of the light curve fit parameters xi and 


V. DISCUSSION 

That the SN la Hubble diagram appears consistent 
with an uniform rate of expansion has been noted ear- 
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Redshift 


FIG. 3. Comparison of the measured distance modulus 
with its expected value for the best ht accelerating uni¬ 
verse (ACDM) and a universe expanding at constant velocity 
(Milne). The error bars include both experimental uncertain¬ 
ties and intrinsic dispersion. The bottom panel shows the 
residuals relative to the Milne model. 


Iiei li 6 | 22 | t 2 g] have confirmed this by a rigorous sta¬ 
tistical analysis, using the JLA catalogue of 740 SN la 
processed by the SALT2 method. We find marginal (i.e. 
< 3 (t) evidence for the widely accepted claim that the 
expansion of the universe is presently accelerating. 


The Bayesian equivalent of this method (a “Bayesian 
Hierarchical Model”) has been presented elsewher^i^. 
We note that a Bayesian consistency tesiP^ has been ap¬ 
plied (albeit using the flawed ‘likelihood’ (equation 12) 
and ‘constrained (equation 13) methods) to deter¬ 
mine the consistency between the SN la data sets ac¬ 
quired with different telescoped. These authors do find 
inconsistencies in the UNION2 catalogue but none in 
JLA. This test had been applied earlier to the UNION2.1 
compilation finding no contamination, but those au- 
thoT^^ fixed the light curve fit ‘nuisance’ parameters, so 
their result is inconclusive. Including a ‘mass step’ cor¬ 


rection for the host galaxies of SN led has little effect. 

While our gaussian model (|^ is not perfect, it appears 
to be an adequate first step towards understanding SN 
la standardisation. One might be concerned that various 
selection effects (e.g. Malmquist bias) affect the data. 
However to address this adequately is beyond the scope 
of this work. We are concerned here solely with per¬ 
forming the statistical analysis in an unbiased manner in 
order to highlight the different conclusion from previous 
analysed^ of the same data. 

We wish to emphasise that whether the expansion rate 
is accelerating or not is a kinematic test and it is sim¬ 
ply for ease of comparison with previous results that we 
choose to show the impact of doing the correct statistical 
analysis in the usual ACDM framework. In particular the 
‘Milne model’ should not be taken literally to mean an 
empty universe since the deceleration due to gravity can 
in principle be countered e.g. by bulk viscosity associated 
with the formation of structure, resulting in expansion at 
approximately constant velocity even in an universe con¬ 
taining matter but no dark energji^. Such a cosmology 
is not prima facie in conflict with observations of the an¬ 
gular scale of fluctuations in the cosmic microwave back¬ 
ground or of baryonic acoustic oscillations, although this 
does require further investigation. In any case, both of 
these are geometric rather than dynamical measures and 
do not provide compelling direct evidence for a cosmo¬ 
logical constant — rather its value is inferred from the 
assumed ‘cosmic sum rule’: Da = 1 ~ Dm + D/j. This 
would be altered if additional terms due to the back re¬ 
action of inhomogeneities are included in the Friedmann 
equation^. 

The CODEX experiment on the European Extremely 
Large Telescope aims to measure the ‘redshift drift’ over 
a 10-15 year period to determine whether the expansion 
rate is really accelerating^. 



-2 0 2 
Pull 


FIG. 4. Distribution of pulls (141 for the best-fit model, com¬ 
pared to a normal distribution. 
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FIGURE LEGENDS 


Fig.l: Distribution of the stretch and colour correction 
parameters in the JLA sampl^^, with gaussians super¬ 
imposed. 

Fig.2: Contour plot of the profile likelihood in the Dm— 
Da plane. We show 1, 2 and 3a contours, regarding all 
other parameters as nuisance parameters, as red dashed 
lines, while the blue lines are 1 and 2cr contours from 
the 10-dimensional parameter space projected on to this 
plane. 



Fig.3: Comparison of the measured distance modulus 
with its expected value for the best ht accelerating uni¬ 
verse (ACDM) and a universe expanding at constant ve¬ 
locity (Milne). The error bars include both experimental 
uncertainties and intrinsic dispersion. The bottom panel 
shows the residuals relative to the Milne model. 


Fig.4: Distribution of pulls (14) for the best-fit model 
compared to a normal distribution. 


Fig.5: The distribution of the likelihood ratio from 
Monte Carlo, with a distribution with 10 d.o.f. super¬ 
imposed. 


METHODS: CONFIDENCE ELLIPSOIDS 


The confidence ellipsoid is the collection of points x = 
{Dm, Da, a, a:o,cr2^,/3, Co, cr^^. Mo, which obey 


[x — xm'L¥]^[x — a;MLE]'’" < (15) 


where is a symmetric matrix and ccmle is the MLE. 
The enclosed volume is a conhdence region with cover¬ 
age probability corresponding with high precision to the 
value obtained from Equation (10). The eigenvectors 
of !F are then the principal axes of the ellipsoid, and 
the eigenvalues are the inverse squares of the lengths of 
the principal axes. We approximate this matrix with 
the sample covariance from the MC of section |IV| as 
F = cov(x, x)~^. 


FIG. 5. The distribution of the likelihood ratio from Monte 
Carlo, with a distribution with 10 d.o.f. superimposed. 

To make reading the matrix of eigenvectors easier, we 
round all numbers to 0.1. Thus, we get the following 
approximate eigenvectors of F, in columns 



/ 

0.5 

0 

0 

0.8 

0.1 

1 

o 

to 

0 

0 

0 

0\ 



0.8 

0 

0 

-0.5 

-0.1 

0.2 

0 

0 

0 

0 

a 


0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

Xq 


0 

0 

0 

-0.1 

1 

0 

0 

0 

0 

0 



0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

p 


0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

Co 


0 

0 

0 

0 

0 

0.1 

0 

1 

0 

0 

^co 


0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

Mo 


-0.1 

0 

0 

0.3 

0.1 

1 

0 

0.1 

0 

0 


V 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0/ 


with respective lengths of semi-axes 


10"^{172, 85.1,49.8,43.9,38.1, 
9.89,5.93,4.24,1.01,0.304} 
We also list the rounded correlation matrix, 

/Dm 
0.9 Da 


0 

0 

a 





0 

0 

0 

Xo 




0 

0 

-0.1 

0 




0 

0 

0 

0 

0 



0.1 

-0.1 

0 

0 

0 

0 

Co 

0 

0 

0 

0 

0 

-0.3 

0 

-0.2 

-0.6 

0 

0 

0 

0.1 

0.2 

0 

0 

-0.1 

0 

0 

-0.3 

0 


Co 

0 


Mo 

0 


(17) 


’Mg/ 
(18) 

We see that the only pronounced correlations are between 
Dm, Da and Mg. This is also apparent from Table |lj 


GODE AVAILABILITY 

The code and data used in the analysis are available 
at: http://dx.doi.org/10.5281/zenodo.34487 
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TABLE I. Maximum likelihood parameters under specific constraints (in boldface). 


(-21ogr„,ax = -214.97) 


Constraint 

-21og/l/£max 


SIa 

a 

*1,0 

^^1.0 

P 

Co 

CCQ 

Mo 

O'Mg 

None (best fit) 

0 

0.341 

0.569 

0.134 

0.038 

0.932 

3.059 

-0.016 

0.071 

-19.052 

0.108 

Flat geometry 

0.147 

0.376 

0.624 

0.135 

0.039 

0.932 

3.060 

-0.016 

0.071 

-19.055 

0.108 

Empty universe 

11.9 

0.000 

0.000 

0.133 

0.034 

0.932 

3.051 

-0.015 

0.071 

-19.014 

0.109 

Non-accelerating 

11.0 

0.068 

0.034 

0.132 

0.033 

0.931 

3.045 

-0.013 

0.071 

-19.006 

0.109 

Matter-less universe 

10.4 

0.000 

0.094 

0.134 

0.036 

0.932 

3.059 

-0.017 

0.071 

-19.032 

0.109 

Einstein-deSitter 

221.97 

1.000 

0.000 

0.123 

0.014 

0.927 

3.039 

0.009 

0.072 

-18.839 

0.125 



